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ABSTRACT 


The local temperature field of a downgoing slab is 
investigated by using a finite-difference numerical 
approach. A numerical model which simulates the downgoing 
Slab is used to. study the effects of various dip angles, 
different rates of subduction, heat sources and rising 
Material from the upper surface of the slap. The model 
assumes a simple descent mechanism near the surface and this 
mechanism is discussed in terms of the associated earthquake 
field and the interactions of the material there. The rate 
of subduction and amount of shear-strain heating along the 
upper surface of the slab are important factors in the 
determination of the thermal regime. When melting occurs, 
rising material from the top of the slab produces high heat 
flow values at the surface of the earth on the continental 
side of the oceanic trench. The gravity effects or the 
downgoing lithosphere are also investigated. The results 
indicate that the presence of rising melt wili mask the 
gravity effect of the cold sinking slab at low subduction 


velocities (i.e. 1 cm/year). 


The behaviour of time-varying electromagnetic fields 
corresponding to various temperature distributions 
associated with a subducting lithospheric slab has been 
studied using a two-dimensional numerical model. The resuits— 
show that the movement of melted material from the top of 
the slab considerably affects the spatiai behaviour ot tne 
electromagnetic field componerts at the surface, 
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particularly that of Hz and H,/Hy. A lateral conductivity 
contrast at the surface (i.e. a sea-land interface) 
completely dominates the behaviour of the field components 
there, and the effect of any subsurface temperature 


variations is negligible. 
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CHAPTER 1 INTRODUCTION 


Associated with the creation of new crustal material at 
ocean ridges and the movement of lithospheric plates is the 
downwarping and descent of plates under island arc regions. 
The concept of ocean ridges and subduction POnee is 
incorporated in the. large-scale kinematic model of the 
earth's Tucan hee known as plate tectonics (Le Pichon et 
al,1973; Isacks et al, 1968). In recent years interest has 
developed in the investigation of the descent of 
lithospheric plates in subduction zones and the island arc 
regions associated with them (Oxburgh and Turcotte, 1968; 
McKenzie and Sclater, 1968; Mckenzie, 1969; Oxburgh and 
Turcotte, 1970; Minear and Toksoz, 1970; and others). Such 
investigations seek to develop a model of the descending 
lithosphere which is consistent with various geophysical and 


geologic data available. 


One of the most important geophysical observations 
taken in recent years has been surface heat flow in relation 
to subduction zones. One of the first heat flow naps 
constructed on a world-scale was that compiled by tLee and 
Uyeda (1965). It included heat flow measurements from ocean 
trench and island arc regions. The ocean trenches were found 
to have heat flow values averaging 0.99+0.61 HOF SUS 
(H.F.U.=microcalories/cm@-sec), substantially lower than the 
surface mean heat flow of 1.5+10% H.F.U. Island arc regions 
had widely scattered heat flow values with a mean of 
1.36+0.54 s.d.e(Sstandard deviation), H.F.U., although a large 
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number of the measurements were above the surface mean heat 
flow. A more up-to-date heat flow map (Von Herzen and Lee, 
1969) indicated abnormally high heat flow values associated 
with island arc regions. Horai and Uyeda (1969) found 
anomalously high heat flow values associated with volcanic 
areas that were related to island arcs. Most heat flow 
Observations taken to date have been found in the Sea of 
Japan where heat flow values average 2.01+.38 s.d., H.F.U. 
(Sugimura and Uyeda, 1973). Since island arc regions in 
general are associated with volcanic areas, it is generally 
believed that most other island arcs display a corresponding 


high heat flow pattern. 


In the past several years Many studies have 
investigated possible causes of anomalous heat flow values 
related to subduction zones as well as thermai regimes 
associated with the descent of the earth's lithosphere in 
these regions. Oxburgh and Turcotte (1968) determined that 
frictional heating may play an important role in generating 
high heat flow in isiand arc regions. According to their 
analytical calculations, production of magmas occurs between 
100 km and 200 km depth, where temperatures are high enough 
on the descending lithosphere*s fault zone to produce 
Significant partial melt. The partial melt was the 
hypothesized to rise to the earth's surface creating an 
anomalous heat flow region. McKenzie and Sclater (1968) 
investigated high heat flow in the northwest Pacific and 


determined that possible sources of heat included volcanic 
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intrusion into the crust and dissipative shear heating in 
the mantle due to the descending lithosphere's motion. 
Caiculations demonstrated that shear heating near the 
earth's surface or an increase in thermal conductivity due 
to mass transport of magma could explain the high heat flow, 


although the initial temperatures needed were unrealistic. 


McKenzie (1969) obtained a temperature distribution in 
the lithosphere thrust beneath an island arc region. Flow 
and stress heating in the mantle maintained a high heat flow 
anomaly in the island arc. The analysis assumed the 
temperature gradient below the lithosphere followed the 
adiabatic gradient and the descending lithosphere acted as a 
thermal boundary layer supporting temperature and density 
gradients larger than those in the surrounding mantle. fhe 
heat flow anomaly in the ‘real world’ situation could then 


be explained if there existed a large shallow stress field 


to produce Shearing near the earth's surface. The 
calculations included a large number of linearized 
approximations. 

Oxburgh and Turcotte (1970) assumed frictional 


dissipation of heat occurred along’the fault zone near the 
upper surface of the descending slab and established a line 
heat source along this zone with source strength as a 
function of depth. They found that most of the heat was 
trapsferred to the descending lithosphere and the heat 
source was unable to generate enough heat to correspond to 


high heat flow in the island arc region. Minear and Toksoz 
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(1970) and Toksoz et al (1971) used a finite-difference 
solution of the conservation of thermal energy equation to 
determine the effects of Spreading rate, adiabatic 
compression, radioactivity, phase changes and strain heating 
on the thermal regime of the downgoing lithospheric’ slab. 
They were unable to generate the heating required to produce 
a heat flow anomaly in the island arc region. Hasebe et al 
(1970) considered two-dimensional time dependent heat 
conduction equations with large heat generation along the 
deep-earthquake fault zone and effective heat transfer due 
to rising magma in the upper mantle above this zone. It was 
determined that heat generation four times aS high as the 
mean heat flow at the upper surface of the lithosphere and 
heat conduction ten times as high as normal phonon 
conduction were required to explain anomalously high heat 
flow in the island arc region. The values used were 
justified by invoking viscous heat generation due to shear 


stress and penetrative convection of magma. 


Griggs (1972) established a thermal model of the 
descending lithosphere which included gravity effects and 
possible earthquake mechanisms related to slab motion. 
Richter (1973) proposed fluid dynamic models of the 
lithosphere-asthenosphere system which included generalized 
mantle convection. He obtained flow patterns within the 
mantle for idealized models assuming the mantle acted as a 
Newtonian fluid. Turcotte and Schubert (1973) determined 


temperatures around the descending lithosphere using 
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constant shear stress on the fault zone as well as stress 
weakly dependent upon depth. The analytic model assumed 
partial melting occurred along the fault zone below the 
point where active volcanoes lie. The temperatures and shear 
stresses agreed with the models proposed by Toksoz et al 


(1971). 


Andrews and Sleep (1974) used numerical modelling of 
flow induced by the downgoing slab to explain high heat flow 
in marginal basins and island arc volcanism. The induced 
flow was two-dimensional, time-dependent and included 
viscosity dependent upon temperature and pressure. The 
highest stress occurred at the base of the slab interface 
creating a narrow channel for separating melt from source 
materials, obtaining volcanic magma. In marginal basins 
tensile yielding caused the plate to thin and raise heat 
flow values. Bird et al (1975) computed thermal regimes for 
continent-continent convergence zones modelled by the 
finite-difference technigue of fToksoz et al (1971). Their 
model of the Zagros mountain region satisfied gravity, 


seismic, heat flow and geologic constraints. 


Spence (1977) established a slabbing, escalator-like 
descent mechanism for the subducting Pacific plate. He 
proposed that partial melt is concentrated in a thin tabular 
zone of mantle above the plate. Partial melt is concentrated 
by episodic pressure reductions due to viscoelastic rebound 
induced by loading of the mantle. Partial melting caused by 


deeper pressure pulses rises to the underthrust plate 
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creating a marginal basin. 


The above studies indicate that there is interest in 
the interaction between the descending slab material and the 
Surrounding medium as well as the form of the downwarping 
mechanism near the surface, though very 1Jlittle has been 
published on this latter point. In the work of Minear and 
Toksoz (1970) as well as Toksoz et al (1971) the thermal 
regime of a 45° dipping slab was investigated. It also 
appears that in most discussions concerning theoretical 
calculations for downgoing slabs, a 45° dip angle has been 
assumed. It is known that the earthquake zones (fault 
zones), which are taken to define the upper boundary of the 
downgoing slab, dip at various angles depending upon depth 
and the region being considered, not necessarily at 45° (see 
Sykes, 1966; Benioff, 1954). The dip angles may vary from 
approximately 15° in the Zagros mountain region and 229° for 
the near surface portion of the slab near the west coast of 
South America to as much as 75° in the deeper Bonin-Honshu 
region. In this regard Bird et al (1975) have recently 
considered slabs with shallower slopes (30° and 159%) and, as 
has been mentioned, have applied their method to the Zagros 


mountain region. 


Tt is also apparent from the profiles of earthquake 
epicenters that the dip angle of many zones changes from a 
smaller angle near the surface to a greater angle at depth 
or that zones may flatten out at greater depth. If we 


identify these zones of earthquake activity with the upper 
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boundary of the slab, the implication is that the moving 
lithospheric plate must ‘bend! at depth as well as initially 


near the earth's surface. 


Several conductivity structures which are postulated to 
account for various regional geomagnetic field anomalies are 
attributed to variations in the temperature below the 
earth's surface (Uyeda and Rikitake, 1970; Law and 
Riddihough, 1971; Gough, 1973, Garland, 1975). Variations in 
the subsurface temperature may be either due to large-scale 
tectonic activity in the crust and mantle regions or they 
may be associated with the subduction of lithospheric plates 
into the upper mantle. Examples of geomagnetic anomalies 
that are attributed to tectonic activity include those 
observed at Alert (Rikitake and Whitham, 1964; Niblett et 
al, 1969), Northern Germany (Schmucker, 1959, Knodel, 1968) 
and Texas (Schmucker, 1964). Rikitake and Whitham (1964) 
suggested that certain features of the Alert anomaly could 
be explained by the rise of the 1400-1500° isotherms to 
within 25 km of the earth's surface. Schmucker (1964) 
attributed the Texas conductivity anomaly to the rise of the 
high temperature isotherms in ‘the upper mantle. This 
interpretation by Schmucker has been supported to a large 
extent by the heat flow observations of Warren et al (1969). 
Comparing their heat flow observations with the geomagnetic 
observations of Schmucker (1964), Warren et al (1969) found 
that regions where Schmucker observed large anomalies in the 


horizontal magnetic component and small anomalies in the 
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vertical magnetic component, had high heat flow near the 


ground surface. 


Two well known examples of geomagnetic anomalies 
observed in sere eon zones are the Japanese anomaly and 
the Andean anomaly of Peru. It has been suggested (Gough, 
1973, Garland, 1975) that both of these anomalies are 
related not only to the downward movement of the cold plate, 
but also the upward movement of hot material produced by the 
frictional heating along the surface of the descending 


lithosphere. 


Although an exact relationship between geomagnetic 
field anomalies and heat flow anomalies is yet to be 
developed, it is believed that the anomalies observed in the 
horizontal magnetic components are more closely correlated 
with heat flow variations than anomalies in the vertical 
magnetic component (Rikitake, 1966). Large values in the 
horizontal components are generally associated with an 


increased heat flow. 


Past studies indicate that the descending slab remains 
cooler than the surrounding mantle to depths of 600 km to 
700 km. The slab experiences shear-strain heating as it 
descends, mainly along the upper surface, although the 
amount of shear-strain heating and its effect remain an area 
of controversy. An even less known quantity is the amount of 
heating provided by phase changes to the interior of the 


Slab as it descends. Evidence for the mechanism of slab 


edt ius0eq wolt t53i- dpsed bod a | sid 
d a c oan . “ae Vw 


ee me 
avilsnons voisjeapamoer +6. “ polqmste: - siied: 
bus yiamomes seodngst add pate eereyws 
sApyod} Sideceens n5d Spt it O24. do 8 


% 
he J 
tes > JOU te 


515 seifibmonts seerd to dtod te dt 
sateiqg tifoo Sd +0. ects} Soiaanee a you 

; : het aT} : ; 
Sa3 ya bsbuborg fpti9st6m Joni To $203 VOM » bre iu dt one | qa 
: mth: re ¥ a 


bahbnss ab odt. to saaaree sat whee - patise 


v 


viseapsaosp nsewied gidenottéiaz $28x2 ae 


st ot -toy ai esiismons wolk test bas BEC 
. | . 7 7 cay 
dy mi bsyisedo astismoas odF sland pevei ied Saath 


i] 
a 
Wye 


beisis1z1059 ylseclo s10M 9th 23 sea0quoo abt 6 


is 


LSoitisv edd ac estismons o5a3 aiveueis cau 
ce ohm bah eR 


od3 af eopisvy sprsd _ 4: oa8h" eaeazaeay aie 


! eu | 1) 


= a : i" ay de iy bee 
as s«dsiw bstsioo0ees ‘yllessaap S26 atAgQgednod 
a a ‘ y A / te ‘ 


4 


\ a 


4 5 


antene7  dsle pit birsde ish. eda sends othobbua eotbuae yeeo 
ot aX Oa 10, aisdas o+ of tacn enibavorawe, ‘ods Mine 
$i as paid sod abbites1kade z nou kagixe ‘dsie od, 


ons Aguodsis B0e2 Ne teqqu eas -phols Hsiss 
6915 Az niga da toatte ey bas Pdi toon histie~rapia, re at 
pie . ) 


to’ tauons od3 a Ytisavsp. awOnd | zee gave) ‘oh. 


} M vad 


947 io Toksodas edt oF syndy. eri 


dois to tien ous 208 coments 


: ed “ A | j Sve ‘a { f ura Ngee a 
n : ay BY 7 ites - 


44 
om 
ae 
wits s4 <> ft 


descent indicates that the slab experiences normal faulting 
as it begins to subduct followed by thrust faulting implying 
the slab acts as a distinct mechanically separate unit from 
the surrounding mantle. The effect that density and 
temperature variations have on the associated gravity and 
electromagnetic fields above descending slabs is still 
uncertain. Current studies also indicate that the descending 


Slab does not maintain a constant dip angle as it descends. 


Considering these aspects of the evidence of slab 
motion, it is of interest to develop a model which can 
provide for a range of dip angles as well as for the more 
general case of a change in dip angle at depth. This gives a 
flexible model which may be used to determine variations in 
geophysical parameters at depth and near the earth's 
surface. It is then possible to investigate the lower 
limiting case of conditions required to produce melting on 
the upper surface of the downgoing Slab. Several theoretical 
models for the thermal evolution of a downgoing slab are 
considered, including varying dip angle, inclusion of 
various heat sources at depth, and rising material from the 
top of the slab. This allows for the development of a model 
which predicts the possible loci of melting at depth and 


variations of heat flow at the surface. 


The properties “Oe” The earth's interior canbe bet cer 
understood if data can be coliected by various geophysical 
methods (i.e. heat flow, gravity, seismic, and 


magnetotelluric measurements) and correlated with each 
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Other. In recent years, studies involving such detailed 
correlation are becoming more common (Garland, ASTS) « 
Nevertheless, geophysical surveys involving the simultaneous 


use of more than one method may not always be practicable. 


One relatively simple means of undertaking a 
correlation study between various methods is through use of 
computer modelling. A numerical model is established to 
compare the thermal properties of a descending slab and the 
associated gravity effects. A comparison is also undertaken 
to determine the relationship between heat flow anomalies 
and electromagnetic field anomalies for a two-dimensional 
model as well as the electromagnetic effects of temperature 
variations within the descending Slab. Electrical 
conductivity models are derived from temperature 


7 


distributions for various thermal models. 
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CHAPTER 2 THE NUMERICAL MODEL 


In the numerical model to be presented, the dynamics of 
lithosphere motion are assumed and perturbations in various 
geophysical quantities are calculated on the basis of the 
dynamics. This numerical method is based on that originally 
developed by Minear and Toksoz (1970). The region of 
interest is divided into a two-dimensional grid with initial 
conditions specified at each grid point. Computation of the 
temperature field proceeds by using the thermal energy 
balance equation and the Peaceman and Rachtord alternatiag- 
direction implicit finite-difference numerical technique 
(Peaceman and Rachford, 1955). The variations in 
temperatures result in density variations which represent a 
variation in the gravity over the region of interest. This 


variation is calculated and plotted. 


In the model presented by Minear and Toksoz (1970) 
novement of lithospheric material into the mantle is 
accomplished by translation of temperatures ahead otf the 
slab (see Minear and Toksoz (1970), Fig. 4). The material, 
as it is being subducted, is assumed to flow ‘around the 
bend* in circular arcs. This method requires a linear 
interpolation procedure to determine the temperature values 
on the face of the slab as the slab bends downwards. 
Material assumed to be at a grid point may not stop at a 
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corresponding grid point on its ‘journey into the nantle 
over the period of time being considered, hence, the _ need 


for linear interpolation between grid points. 


In the method presented here, the downward slab 
movement 1S carried out by translating the temperature 
values diagonally through the mesh as illustrated in Fig. 
1(a). This is advantageous in that no interpolation is 
necessary. Also, different dip angles for the Slab‘s descent 
may be readily accommodated by varying the relative lengths 
of Ax and Az as shown in Figs. 1(b) and i(c). This method 
can also conveniently provide for a change in the dip of the 


Slab as illustrated in Fig. 1(c). 


Comparison of the present model with that of Minear and 
Toksoz leads to considerations of the mechanism by which the 
slab begins to subduct and the nature of the interaction 
between the moving oceanic plate and the stable plate under 
which the slab is descending. Whereas the Minear and Toksoz 
model implies a bending of the slab by flexural slip with 
tensional stress at the top of the bend in the slab and 
compressional stress at the bottom of the bend, this model 
implies a downward slipping of the end of the slab together 
with local tension as it starts to subduct, essentially due 
to a normal faulting mechanism. Consider an initial steady- 
state lithosphere, becoming gravitationally unstable over 
lower density material: it will sag at first, creating a 
broad trench-like region on the ocean floor, followed by an 


inverted keystone failure at the axis of flexure. Successive 
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normal faults will propagate into the slab, parallel to the 
first, permitting its end to descend into the mantle. The 
analogy can be made with a pack of cards sliding past one 
another. The stress-field at the top of the slab will 
necessarily be tensional subparallel to the surface. 
Calculations done by Lliboutry (1969) indicate that the 
stresses involved in lithospheric sinking allow for a Simple 
Shear mechanism with normal faulting. This model maintains 
the vertical thickness of the lithosphere is unchanged with 
downward sinking of the plate in blocks by shearing along 


vertical planes. 


This type of lithospheric descent is supported by first 
motion studies done by Isacks et al (1968) which indicate 
stress and deformation within a lithospheric plate may occur 
by normal faulting at shallow depths beneath the axis of the 
trench and landward of the underthrust zone. This type of 
normal faulting beneath the seaward slope of the oceanic 
trench is further supported by Katsumata and Sykes (1969) 
with studies done in the Phillipine Sea. Malahoff (1970) has 
concluded that gravity faulting is the dominant type of 
deformation in several trenches in the Pacific based on data 
obtained from seismic reflection profiles. Ludwig et al 
(1966) believes gravity faulting to be dominant in the Japan 
trench and Stauder (1968a, b) reached similar conclusions 
for the Aleutian trench. Lister (1971) obtained on acoustic 
reflection prctile across the Chile trench off Valparaiso 


which showed downbowed reflections and normal faulting 
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indicating the lithospheric slab may . be descending in an 
escalatorlike fashion. In terms of the probable history of 
the slab, the above translational procedure is believed to 
represent a possible mechanism for lithospheric movement in 


subduction zones. 


Calculation of the temperature field is accomplished 
using the Peaceman and Rachford alternating-direction 
implicit finite-difference method (Peaceman and Rachford, 


1955) and the thermal energy balance equation: 


aT oz. ‘ 
Ca ne (1) 


where Cp is the specific heat at constant pressure, Pp is the 
density, T is the temperature in degrees Celsius, K is the 
thermal conductivity and H is the heat generation per unit 
volume. Computation of the thermal regime proceeds as 
follows: temperatures representing lithospheric material are 
translated downwards as described in the previous’ section; 
the thermal energy balance equation is then employed to 
obtain the corresponding temperatures after a time At; the 
lithospheric motion is then resumed and the temperatures 
recalculated; and so on. The vertical translation distance 


for the lithospheric slab is 10 km for each time step. 


Use of the thermal energy balance equation is 
accomplished by expressing the equation in finite-difference 


form and using an algebraic scheme proposed by Thomas 


as ak YakBobest ed yee dske bia 
to Yrosakd sidstorg eds" To" anaae | 
ot hewsiled ai ain hesorg. aes 
nk t7enevoa pizedqeodsil rox" ueinsdos 


; 


oistigso2 bheH Pere 


. ay 


at, 
bywintigaooss si blot siutbIsqme3 ie 
aoitooxLb-pgisenistis bro? dosh bats. a 
»bsoOTADSHA Dh éeiosesey fiditea eono1eX1ib~9F, 


4 ; 
Lie 
[ ' 7: 


:HoitsEps 991 51.50 yexsas Lanzed 


st : \ . We 
zi ‘ . ae lao ¥ 
( r) H ot (TV Sh = ed fe i a a i i 


oAt AL Q wIhezere te she BOS: (th yesd DERROSge 
$a at A ,eutelod 20oxpep nt ois $niegasd oat: 
sayy 


fino <teg cade anes $s5d° onF ar a bas yiive: 


25 ia snipes Jsmarsat oe fia is abe 


mots9sa pach ly ant ae bodiaseebh 38 Maa 
oF become, aedd 25 HOLIEHpS , soisted. (wasn Lem doa aia 


edt SA ML 6 s09%6 29m $3193 firs “pnibaogeo 7x0 ony per 
eoLut 679g mes! odd” bas bemvess nad pi ‘PeRtOR bis au 


90% bz cb HoitsieasTs PASt aise ad® ino oe bas ong 


<qod2 omit sobs soY mA OF ai ini » ohvodaeedgn ait a ie 


) i B a" ; z 


ef uot supe ~ sonsied yex909  Kemxaitd - prem cailil a ae 
- bebe i i 


- 


e007 922 LbSS Fink? At S017 6tpe oa ; 


) 


nik 


sa2029%9 “yd: bawaiks 
lead ‘ "oi Pe a 
oh 


eemod't a bézogoiy “auedoe dkeqdapte a5: vaien ‘ ieiiaanee: As 2 


16 


(Peaceman and Rachford, 1955) and modified by Minear and 
Toksoz (1970) to obtain temperature solutions for each time 
step. The scheme requires setting up the derivatives to 
obtain temperature solutions solved implicitly in the x- and 
z-directions in alternate steps, resulting in the final 
solution. A comparison between the solutions obtained from 
the finite-difference method and an exact solution for _ the 
thermal energy balance equation indicates that the solutions 


agree to within a few percent (Toksoz et al, 1971). 


Stability criteria established by Peaceman and Rachford 
(1955) indicate. that the finite-difference solution is 
unconditionally stable for the thermal energy balance 
equation for the case of constant thermal conductivity. In 
other words, the scheme converges to the correct oats on 
for all grid sizes and time steps. In the case of thermal 
conductivity being temperature dependent, the stability 
conditions have been established by Toksoz et al (1971) and 
Minear and Toksoz (1971). The stability conditions require 
that the spatial derivatives in thermal conductivity be 
restricted in magnitude to allow convergence of the scheme 
to the proper solution. In all examples presented, the 


convergence requirements have been met. 


The boundary conditions are similar to those employed 
by Minear and Toksoz (1970). Heat flux into the bottom of 


the grid was taken as constant to maintain steady-state 
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conditions in the mantle. The surface of the earth is 
maintained at 0°C and it is assumed that there is no heat 
flux across the side boundaries of the region. This latter 
condition implies that the edges of the region of interest 
are not affected by the downgoing slab. The density 
distribution used in the calculations closely follows the 


Bullen A density curve and iS Shown in Fig. 2. 


The initial vertical temperature distribution was 
obtained using the relationship given by Mercier and Carter 


(1975) for low-temperature oceanic pyroxene geotherms: 


T = 4,34(P+8.6) - 11840/(P+8.6) + 1340 (2) 


where T is the temperature in degrees Celsius and P is the 
pressure in kilobars (where P follows the normal hydrostatic 
gradient). Small variations in the initial temperature 
distribution are unimportant Since it is of interest in this 
model to study perturbations in the thermal regime. Also, 
there is much uncertainty in the exact mantle geotherm as 
illustrated by the variety of geotherms published (Mercier 


and Carter, 1975). 


The specific heat at constant pressure is taken as 
constant and given the value Cp = 1.3x107ergs/gm°C (Toksoz 
et al, 1971). Fcr comparisons with the Minear and Toksoz 
method the thermal conductivity used is given by MacDonald 
(1959), which includes constant lattice conductivity plus a 
radiative transfer term which is dependent upon temperature. 


This thermal conductivity gives unrealistic values for 
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19 
temperatures above 1000°C and will therefore affect the 
temperature distribution within the descending slab 
(Turcotte and Oxburgh, 1972). Furthermore, experimental 
evidence indicates the radiative opacity of riianleaane 
increases with temperature (Aronson et ak. ¢ 1970). 
Therefore, for most models a constant conductivity is used 


with a value of K= .01 cal/ycm-sec°C (Hasebe et al, 1970). 


2-4 Heat Sources 


The inclusion of heat source coefficients is 
incorporated through the heat production term in the thermal 
energy balance equation. Heat sources which have been 
included are those due to adiabatic compression and shear- 
strain heating. Heating due to radiogenic sources has not 
been included since it is expected the effect will be 
Minimal in the time periods considered. This has been 
supported by calculations done by Minear and Toksoz (1970) 
and Toksoz et al (1971). Heating due to phase changes within 
the descending slab is of no known importance for depths 
less than 300 km and has not been considered because of the 
uncertainty in the contribution this term has in the 
temperature field computation. Calculations done by Oxburgh 
and Turcotte (1970) indicate the effect of phase changes is 


insignificant on the thermal regime. 


The heating of the lithosphere as it descends into the 
mantle due to compression is calculated by the rate of 


energy release at depth h (Toksoz et al, 1971; Hanks and 
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Whitcomb, 1971): 


dQ Le ‘Oo oT = TV (3) 
dt h if ©,° ot CP Je, ae Zz 


where p is the density, Cp is the specific heat at constant 
pressure, Vz is the verticai velocity of the slab, g is 
gravity and T is the. temperature. The volume coefficient of 


expansion is assumed to have a depth dependence given by: 


a = exp(3.58 - 0.00722) (4) 
with the depth z in kilometers and © in units of 10 ©/°C. 


The inclusion of a shear-strain heating term is 
constrained geologically by the need to create partial melt 
at the upper surface of the downgoing slab (Hatherton and 
Dickinson, 1969). Several shear-strain coefficients were 
employed in order to determine if partial melt was produced 
for the different subduction velocities. Shear heating was 
included along the upper surface of the slab and also along 
the bottom edge of the slab. To estimate the effect of 
shear-strain heating, assuming viscous dissipation: 


2 


es ec Sone (5) 
E=oc€é#= n($7) 
where o is the stress, € is the strain rate, nis the 
viscosity, and the velocity gradient { ss ) «as taken 


perpendicular to slab motion, across the shear zone. If we 
assume the velocity distribution across the shear zone 


varies linearly as: 


v(é) = Yo = ests (6) 
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where a is the slab velocity, T is the width of the shear 
zone and § 1S measured across the shear zone,the rate of 
viscous dissipation is then: 


awe 
7 a mus) (7) 
E (Z| 


According to theoretical caiculations (T. Spanos, private 
communication) the shear zone is approximately 10 km wide. 
Therefore, the viscosity of the shear zone may be estimated 
if the amount of shear heating needed to produce partial 


melt is known. 


The model may be adapted to study the effect of 
generating melt along the upper surface of the downgoing 
Slab and causing it to rise. This simulates rising diapiric 
material when it has reached a stage of partiai melting at 
the top edge of the slab, i.e. when the geotherm has 
intersected the solidus curve for basalt. The solidus for 
basalt was taken from Yoder and Tilley (1962) with a 
gradient of 39C/km. Geological studies indicate that basalt 
magma is generated along the upper surface of the slab and 
this magma undergoes crystal fractionation forming andesitic 
magma near the earth's surface (Hatherton and Dickinson, 
1969; Kuno, 1968). Therefore, simulating upward movement of 


basalt magmas is a good approximation geologically. 


To simulate movement of magma from the upper surface of 


the descending slab temperatures are translated upwards in a 
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Similar manner to movement of lithospheric tateriai 
downwards. The rising material is assumed to have a velocity 
of 1 cMm/year to maintain order of magnitude energy balance 
between heat being produced in the shear zone and the amount 
of heat removed. The amount of heat being removed may be 


estimated by (Hasebe et al, 1970): 
Q = CeveAT (8) 


where C iS the specific heat at constant volume, V is the 
velocity of the rising material and AT is the temperature 
difference between the top and bottom of the magma layer. No 
allowance has been made for heat of fusion as the amount of 
heat required iS about two orders of magnitude iower than 
the amount of heat produced, minimizing the eitfect of heat 


absorption due to heat of fusion. 


2-6 Gravity Calculations 


The gravity effect of a sinking slab is  caiculated 
using a method developed by Dyrelius and Vogel (1971). For a 
two-dimensional rectangular block with upper edge at the 


surface z=0, the gravity effect is given by: 


em acez- sist ll 135 
g = vo{a ln cmeace BD in seer dz (tan Stan 7) U3) 
a b 


where v is the gravitational constant, p is the density, 
x, is the x-coordinate of the heft: edge of the block, x, is 
the x-coordinate of the right edge of the block, ais x-x,, 


bas x-x, and z is the depth to the lower edge of the Dbliock. 
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The above equation also holds for the upper corners of the 


block. When x=x, : 


1 
a [2 In ae = 0 (10) 
and therefore: 
Ad cinaraee [-» In Be 2ztan + | Gig) 
ENG So ae 
AG,-g = VP (= In oe + 2ztan + (#2) 


If the gravity effect is considered from the jth biock at 
the jth point of observation, the total gravity effect at 
the point of observation Ps is: 
m 
Seri, 


i eae (13) 


Jet 


where there are ‘m* blocks. 


For the calculation of the gravity effect a two- 
dimensional grid was used which incorporated the densities 
used in the thermal regime computation. A system of grid 
cells was established which used densities specified at the 
four grid points enveloping a grid cell. The gravity 
contribution trom each grid cell was then added to obtain 
the gravity effect at selected observation points at the 


surface. 


Variations in density that are accounted for in this 
calculation are due to temperature fluctuations in the 


lithosphere and surrounding mantle. The changes in density 
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include the effect of the volume coefficient of expansion 
for temperature above and below the normal mantle geotherm. 
It is expected that the cold sinking slab will hen teo oe a 
higher density than the surrounding material and partially 
molten material will be less dense then the surrounding 
medium. No allowance has been made for density variations 


due to phase changes and chemical inhomogeneities. 
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CHAPTER 3 THERMAL REGIMES FOR DOWNGOING SLABS 


n of the Diagonal Translation method and the 


cers ——eee Sa a oe ee ee ee SoS Soe eaemess Soe ae 


In order to ee ee ae the diagonal translation method 
with the Minear and Toksoz interpolation procedure, a slab 
with a dip angle of 45° was modelled using both techniques. 
A slab 80 km thick superimposed on a grid of 151 x 151 
points Ws used, in which the mesh size is taken as 5 km x 5 
km. The slab moves with a vertical velocity of 1.4 cm/year, 
and so the velocity of the slab through the surrounding 


material is 2.0 cm/year. 


The two models are compared in Figs. 3-5. The left 
sides of the figures show the results for the Minear and 
Toksoz method (Figs. 3A, 4A and 5A), while the right sides 
show the effects of the diagonal translation procedure 
(Figs. 3B, 4B and 5B). The temperature fields are shown in 
the lower diagrams with the corresponding heat flow profiles 
shown at the top of che figures. Figure 3 shows the thermal 
regime and surface heat flow at 10.63 million years after 
the start of subduction, while Figs. 4 and 5 are at 21.26 
million years and 35.44 million years, respectively, after 


the beginning of subduction. 


From the figures the effect of the narrower slab in the 
diagonal translation procedure is apparent in that the 


isotherms lie closer together inside such a slab. This 
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implies that the distortion of the normal geothermal field 
by the cool downgoing slab depends on the thickness of the 
Slab, and the perturbation of the thermal regime throughout 
the surrounding material takes place slowly. Also, the 
effect of the vertical lower face of the slab in the 
diagonal translation procedure is evident, causing 
distortion of the isotherms near this region. Although the 
isotherms differ to some extent near the bottom of the slab, 
the temperature field above the slab is not greatly aftected 
except near the surface. Close inspection of Figs 3 ‘shows 
that the isotherms are slightly different near the surface 
close to the Sorae where subduction begins. This is 
particularly true for the 100°C isotherm. This difference is 
also apparent in Figs. 4 and 5. The surface heat flow vaiues 
reflect this near surface difference and it appears that the 
nininun in the heat flow is lower in the diagonal 
translation method. This difference becomes more apparent as 
time goes on (Figs. 4 and 5) and the minimum for the 
diagonal translation method is noticeably lower at 35.44 


million years. 


The reason for the lower minimum becomes apparent when 
the geometry of the two subducting mechanisms is considered. 
In the Minear and Toksoz method where material is assumed to 
flow ‘around the bend*, the actual vertical motion of the 
material is not as pronounced near the surface of the 
lithosphere as in the diagonal translation procedure. The 


‘bending phenomenon allows the slab to traverse a larger 
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horizontal distance in the Minear and fToksoz method, 
producing smaller temperature gradients near the surface of 
the slab and thus allowing the near surface portion of the 
Slab to warm up at a quicker rate. The effect also causes 
the heat flow profiles to differ as time proceeds, 
particularly to the right of the point -where subduction 


begins (as in Figs. 3, 4 and 5). 


with Different Dip Angles 


The model was used to investigate the temperature 
fields and heat flow profiles associated with a shallower 
dip angle and with a slab in which the dip angle changes 
with depth. The thermal conductivity used in these models is 


given by MacDonald (1959). 


3.2.1 Dip Angle of 26.62 


The results for a slab in which the dip angle is 26.69 
are shown in Figs. 6-8. In this model Ax=10 km and Az=5 kn. 
The vertical velocity of the slab is 1.4 cm/year and 
therefore the subduction velocity is 3.2 cm/year. Figures 6- 
8 show the associated thermal fields and surface heat flow 
profiles of 10.63,  #21.26, and 35.44 million years, 
respectively, after the start of subduction. Comparison of 
this shallower angle with the 45° angle of Figs. 3-5 shows 
the isotherms more severely distorted for the 26.6° dipping 
case. This is explained by the greater subduction velocity 


in the shallower dipping case. The 26.6° dipping slab does 
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not have as much time to warm up as the 45° case, resulting 
in isotherms being taken to greater depths. Also, the 
Shallower dipping slab has a larger thickness in the 
direction of subduction, resulting in a thermal boundary 
layer which requires a greater time to be heated to mantle 


temperatures. 


Although the 45° dipping slab in Figs. 3-5 has a fairly 
Symmetric heat flow profile, the 26.6° dipping slab has an 
asymmetric form which becomes more pronounced as subduction 
progresses. Further, near the point where subduction begins, 
the botton on the dip in the heat flow becomes flatter as 
time progresses (Figs. 6-8). This result 1S apparent from 
the form of the downwarping mechanism and the rate of 
subduction. Whereas the horizontal slab velocity in the 45° 
case is 1.4 cm/year, the horizontal velocity in the 26.69 
case is 2.8 cm/year resulting in a broad minimum in the heat 


flow pattern. 


To investigate the temperature fields of slabs in which 
the dip angle changes with depth, two models were run. The 
first model is that of a slab which dips initially at 26.69 
and after 10.63 million years the dip angle changes to 45°, 
The second model also initially dips at 26.6°, but after 
10.63 million years the dip angle changes to 56.39. In both 
models the downgoing velocity of the material in the 260.69 


portion is 3.2 cm/year. In the lower portions of the slabs 
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the subduction velocity is 4.0 cm/year for the 45° slab and 


5.1 cm/year for the 56.3° slab. 


Figures 9 and 10 give results for these two models at 
21.26 million years after the beginning of subduction. By 
comparing Figs. 9 and 10, it is seen that the slab has 
penetrated deeper into the mantle in Fig. 10 for the same 
time period, resulting in the isotherms being pulied down 
further into the mantle with the steeper dipping slab. The 
lowering of the isotherms can also be interpreted as due to 
the greater subduction velocity in the steeper dipping slab. 
Near the point where subduction begins the heat flow 
profiles and thermal regimes are Similar. This implies that 
the perturbations created at greater depths do not affect 
the thermal regimes near the surface. When compared with the 
slab with constant 26.6° dip angle (Fig. 7) there is littie 
difference between the surface heat flow profiles, although 
the temperatures near the slabs display marked differences 
deeper in the mantle. These results indicate that the 
perturbations in the thermal fields at depth do not have 
sufficient time to be propagated to the surface in the time 
periods considered, giving little information regarding slab 
processes related to heat conduction and slab velocity 


changes at depth. 


In the examples presented here heat conduction along 


with the inclusion of heat source coefficients for shear- 
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Fig. 9. A downgoing slab with dip angles 26.6° 
and 45,0°. 
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Strain heating and adiabatic compression are considered. The 
model was used to compare thermal regimes and heat flow 
profiles associated with slabs dipping at 26.69 and 45° for 
two vertical subduction velocities, 0.7 cm/year and 5.6 


cm/yeac. 


3.3.1 Subduction Rate of 0.7 cm/year 


The results for a downgoing slab with dip angle 45° and 
a Slab with dip angle 26.6° are shown in Figs. 11-13. Figure 
11 illustrates the thermal regime and surface heat flow 
profile associated with a slab dipping at 45° witha 
vertical velocity of 0.7 cm/year. The resulting subduction 
velocity is 1.0 cm/year. Adiabatic compression is included 
as a heat source within the slab together with a  shear- 
strain heating coefficient of 1.6x10-*ergs/cm3-sec along the 
top edge of the slab and 1.6x10 -Sergs/cm3-sec along the 
bottom edge of the slab. The order of magnitude lower shear 
heating term along the bottom edge of the slab is due to an 
assumed smaller viscosity in the low velocity zone in the 
Mantle (Schubert et al, 1976). The shear zones along the top 
and bottom ‘edges of the slab are assumed to have a 10 km 
vertical thickness. Together with the assumed values of 
shear heating, this implies an average viscosity in the 
shear zone along the upper edge of the slab of approximately 
1023 poise. This agrees within an order of magnitude with 
the average viscosity in the shear zone if the descending 


lithosphere is assumed to have a viscosity of between 102% 
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The temperature regime and surface heat flow profile 
for a descending slab with a dip angle of 45° at 
70.88 million years after the start of subduction. 
Adiabatic compression and shear heating are included 
as heat sources. The slab velocity is 1.0 cm/year. 
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Fig. 12. Temperature-depth profiles along the upper surface of 
the slab for the models of Fig. 11 (curve A) and Fig. 
13 (curve B). 
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poise to 1025 poise (Walcott, 1970; De Braemacker, 1977) and 
the upper mantle has a viscosity of approximately 102! poise 
(Schubert et al, 1976). Results for a total subduction 
period of 70.88 million years are illustrated in Figs. 11- 
13. Figure 12 (curve A) gives a temperature-depth profile 
along the upper surface of the slab. This may be compared 
with Fig. 12 (curve C) which represents the soldidus for 


basalt (Yoder and Tilley, 1962). 


From Fig. 11 it is seen that the interior of the slab 
remains cooler than the surrounding mantle to depths greater 
than 500 km jmplying that the slab does not reach thermal 
equilibrium in the time period considered. This may be 
justified by the fact that earthquakes occur to depths of 
700 km within the downgoing slab {Benioff, 1954; Isacks et 
al, 1968] and possibie sources of earthquakes may be due to 
a colder brittle slab descending into the mantle (Griggs, 
1972; Mckenzie, 1969). The shear-strain heating term raises 
the isotherms locally near the top edge of the slab although 
the bottom edge of the slab is unaffected by shear heating. 
The surface heat flow pattern exhibits a distinct low near 
the point of subduction with no visible change in the heat 
flow as a result of shear heating along the upper surface of 
the slab. This result indicates that a thermal pulse 
propagating from the upper surface of the slab will not 
reach the surface in the time period considered. 
Temperatures along the upper surface of the slab intersect 


the melting curve between the depths of 250 km and 350 km 
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implying that melt is generated in this region, 


The thermal regime and surface heat flow profile for a 
Slab dipping at 26.6° are illustrated in Fig. 13. The 
subduction velocity is 1.6 cm/year and the heat sources are 
the same as those of Fig. 11. Total time of subduction is 
70.88 million years. Fig. 12 (curve B) gives the associated 
temperature-depth profile aiong the upper surface of the 
Slab indicating that temperatures along the upper edge of 
the slab reach the partial melt curve from approximately 150 


km to 400 km depth. 


A comparison of the results for the two different dip 
angles indicates that the interior and lower parts of the 
Slab remain cooler at shallower angles of dip due to the 
greater subduction velocity for the shallower angle slab and 
the greater thickness of this slab resulting from the 
descent mechanism. The upper surface of the slab is 
maintained at higher temperatures in the shallower dipping 
slab due to increased shear zone thickness because, in turn, 
of the increased slab thickness at shallower dips. Melting 
then occurs at sShallower depths. The slightly larger 
subduction velocity in the 26.6° slab tends to cause a 
general lowering of the isotherms within the slab region 
although the effect is largely negated along the slab's 
upper surface due to the increased shear zone thickness. 
There is little indication in the surface heat flow of the 
presence of shear heating in either of the slabs, although 


the dip in the heat flow profiles may be more symmetric as a 
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result of shear heating (induced artificially by the depth 


at which shear heating begins in this model: 30 km). 


The effect of generating melt and causing it to rise is 
Shown in Figs. 14-15. For both a 45° dipping slab and 26.6° 
dipping slab containing the same heat sources as shown in 
Figs. 11 and 13, melt was forced to rise for depths 
extending from 250 km to 350 km This simulates rising 
diapiric material after it has reached the basalt solidus 
along the top edge of the slab. The region of rising melt 
extends from the 300 km point to the 400 km point along the 
horizontal axis in Fig. 14 and from the 600 km point to the 
800 km point alone the horizontal axis in Fig. 15. The 
region of rising partial melt is in agreement with studies 
done by Hatherton and Dickinson (1969) on the relationship 
between the depth of the shear zone below volcanoes in 
relation to the distance of the volcanoes from the trench 
axis. A list of a number of volcanoes in relation to trench 


systems is shown in Table 1 (Hatherton and Dickinson, 1969). 


For the 45° dipping slab the result of the melting 
anomaly is to increase the heat flow by approximately 25% at 
300 km from the point of subduction. There is a general 
raising of the geotherms in the melt region due to the 
rising material with a large temperature gradient occurring 
at approximately 30 km depth. Temperatures throughout the 
melt region remain high compared with the surrounding 
mantle. Also, the width of the heat flow anomaly corresponds 


with the width of the melting anomaly resuiting from the 
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Fig. 14. The temperature regime and surface heat flow profile 
for. a,45" dipping slab with rising melt. Heat sources 
are the same as in Fig. ll. 
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Fig. 15. As Fig. 14, 


but with a dip angle of 26.6°. 
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Table 1. Geometry of Active Volcanoes 


am qameqe eas a Temp enarmerae= eee asaemam 


Indonesia No.t h, : km ad, * km 
(21) Marapi 6,1-14 160 320 
(22) Tandikat 6,1-15 130 300 
(23) Galunggung 6,3-14 130 270 
(24) Tjeremai 6,3-17 210 320 
(25) Slamet 6, 13-18 210 ~ 
(26) Dieng 6, 3-20 260 330 
(27) Ungaran 6, 3-23 290 — 
(28) Merapi 6, 3-25 225 _ 
(29) Paluweh 6,4-15 240 _ 
(30) Lewotolo 6,4-23 240 _ 
(31) Lokon-Empung 6,6-10 130 oe 
(32) Dukono 6,8-1 180 _ 

Lesser Antilles 
(33) Mt. Misery 16-3 120 250 
(34) Nevis Peak 16-4 120 250 
(35) Montserrat 16-5 120 250 
(36) Mont Pelee 16-12 130 _ 
(37) Qualibou, St. Lucia 16-14 130 ~ 

New Zealand 
(38) Ohakune-Edgecuombe ~ 140 300 
(39) White Island _ 210 300 
(40) Egmont - 280 380 


Number of volcanoes in Catalogue of Active Volcanoes 


of the World [1951, 


1961]. 


depth to center of Benioff zone below volcano 


distance to volcano from centre of trench, 


present. 
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vertical transport of magma (and therefore heat) towards the 
earth's surface. Figure 15 shows the melting anomaly in the 
26.6° dipping case rising from the same depths as those of 
Fig. 14. Due to the revised grid dimensions for the 26.69 
dipping Slab, the region of partial melting is widened along 
the horizontal axis. Heat flow is also raised by 
approximately 25% over the width of the anomaly. The region 
of high heat flow is in agreement with heat flow 
measurements in the Japan Sea (Vacquier et al, 1966) which 
indicate that the width of the heat flow anomaly in the 
island arc region is dependent on the dip angle of the fault 
zone below the island arc. A corresponding high heat flow 
anomaly over a shallow dipping slab (approximately 30° dip) 
has twice the horizontal width as compared with a steeper 
dipping slab (approximately 45°). This result may indicate 
that the corresponding heat flow anomaly due to partiai 
melting may ke caused by melt created in a specific depth 


interval not related to the dip angle of the slab. 


The results for a downgoing slab with a vertical 
velocity of 5.6 cm/year are Shown in Figs. 16-19. In the 45° 
dipping slab of Fig. 16, the subduction velocity of the slab 
is 8.0 cm/year. This represents an ainper limit for the 
spreading rate of the lithosphere (Le Pichon et al, 1973). 


Total time of subduction is 8.86 million years and heat 
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Fig. 16. The temperature regime and surface heat flow profile 


for a 45° dipping slab with a velocity of 8.0 cm/year. 
Heat transport is by conduction only. 
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conduction is considered as the only process which governs 
the slab's thermodynamics. Figure 16 shows the associated 
temperature field and surface heat flow profile for the 
descending slab. From the figure it is seen that the slab 
remains cooler than the surrounding mantle to depths greater 
than 500 km indicating that heat conduction alone is 
Uy ouieai ode nt to allow slab assimilation into the upper 
mantle. The increased slab velocity as compared with that of 
Figs. 3-5 results in the heat flow profile dipping to near 
zero values at the point of subduction. This suggests’ slab 
dynamics may be indicated by the magnitude of heat flow 
observations in ocean trenches. Lower heat flow values in 
one trench relative to another may indicate that the 
velocity of the slab with lower heat flow may be greater 
than the slab which displays higher heat flow in the trench 
region. It is uncertain whether this effect will be masked 


by regional anomalies of a near surface origin. 


The effects of adiabatic compression and shear-strain 
heating on a 45° dipping slab with a subduction velocity of 
8.0 cm/year are shown in Fig. 17. Shear heating occurs in 
the same areaS as mentioned with regard to Fig. 11 and has a 
magnitude of 6.0x10-* ergs/cm3-sec along the top surface of 
the slab and an order of magnitude lower term along the 
bottom surface of the slab. The implied viscosity in the 
shear zone along the top Surface “of .the slab is 
approximately 1022 poise. The combined effect of adiabatic 


compression and shear heating is to warm the interior of the 
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Fig. 17. As Fig. 16, but with adiabatic compression and 
shear heating included as heat sources. 
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Slab substantially as compared with heat conduction alone. 
Shear heating along the upper surface of the slab raises the 
temperatures above the basalt solidus (Fig. 18, curve C) in 
depths ranging from 175 km to 350 km (Fig. 18, curve A). The 
effect of shear heating is not observed at the earth's 
surface due to the periods of time considered. Figure 18 
(curve B) gives the temperature-depth profile (along the 
upper surface of the slab) for a slab dipping at 26.69%. The 
Shear zone thickness for the 26.69 slab is slightly larger 
than that for the 45° dipping slab in the direction of 
subduction, resulting in more intense heating along the 
Shallower slab. This increased shear zone thickness is a 
result of the descent mechanism as described previously 
(Sec. 3.3.a). The vertical velocity of the 26.6° dipping 
Slab is 5.6 cm/year, resulting in a subduction velocity of 
12.5 cm/year. The effect of a much larger’ subduction 
velocity is negated by shear heating along the slab"s upper 
surface since it is expected a greater subduction velocity 
will tend to maintain cooler temperatures aloag the slabs 
upper surface. The region of partial melting is extended 
from 100 km to approximately 400 km depth in the 26.6° case 


(Fig. 18, curve B). 


The effect of increased shear-strain heating in the 
26.69 and 45° dipping Laws is illustrated in Fig. 19. 
Figure 19 (curve A) gives the temperature-depth profile 
along the upper surface of a 45° dipping slab with shear 


heating along this upper surface amounting to 7.0x10-¢ 
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Fig. 18. Temperature-depth profiles along the upper surface of 
the slab for Fig. 17 (curve A) and for a slab dipping 
at 26.6° with the same heat sources as in Pig .al 7 
(curve B). 
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Fig. 19. As Fig. 18, but with shear heating increased by 15% 
for a 45° dipping slab (curve A) and a 26.6° dipping 
slab (curve B). Curve C is the basalt solidus. 
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ergs/cm3-sec. The corresponding 15% increase in shear 
heating along the upper surface extends the region of 
partial melting from 100 km to 400 km depth, an increase of 
70%. A similar increase is noted for the 26.69 case as shown 
in Fig. 19 (curve B). The shear zone has aé_e shear heating 
coefficient of 7.0x10-* ergs/cm3-sec and is 10 km in 
vertical thickness along the upper surface of the slab. The 
thickness of the region of partial melting is increased by 
20%. The above results indicate that a small increase in the 
Shear heating along the upper surface of the slab produces 
large changes in the thermal fields of the respective slabs. 
Shear heating may control the regions of partial melting in 
the descending slab while being directly controlled itself 
by the amount of stress propagated by the coupling of the 


Slab-mantle interface. 


3.4 Gravity Effect 


Figures 20-21 indicate gravity effects due to sinking 
slabs with various heat sources and rising material from the 
upper surface of the slab. Figure 20 illustrates the gravity 
effect for slabs dipping at 45°. In Fig. 20, curve A gives 
the results for a slab with a vertical velocity of 
taubaveribin of 0.7 cm/year and adiabatic compression included 
with shear heating along the upper surface of the slab (as 
in Fig. 11). The resulting gravity effect indicates a weak 
negative anomaly over the horizontal axis. The effect of the 


cold sinking slab is masked by the partial melt created at 
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Gravity profiles across the horizontal extent of 


subduction for a 45° dipping slab. Curves Ae BS TC 


and D correspond to the models of Figs. 11, 14, 16 
and 17 respectively. 
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the upper surface of the slab. Figure 20 (curve B) 
represents the gravity effect due to a 45° dipping slab with 
partial melt rising from the slab‘'s upper surface (as in 
Fig. 14). The resulting gravity anomaly shows a_ broad 
regional low with a trough of approximately 60 milligals. 
The effect of the rising melt again masks the dynamics of 
the cold descending lithosphere. The anomaly has a long 
wavelength which results from the horizontal extent over 
which rising melt occurs. Figure 20 (curve C) illustrates 
the effect of a 45° dipping slab with subduction velocity of 
8.0 cm/year (as in Fig. 16). The effect of a cold sinking 
Slab warmed only by heat conduction is represented by a 
positive regional anomaly with a peak of approximately 75 
Milligals. The large positive anomaly results from the large 
subduction velocity causing the lithosphere to remain cooler 
than the surrounding mantle to depth of 500 km. In Fig. 20, 
curve D is the gravity effect of an 8.0 cm/year subducting 
Slab warmed by adiabatic compression and shear heating (as 
in Fig. 17). The additional heat sources cause the peak of 
the anomaly to be lowered to about 25 milligals as a result 
of less dense material being formed along the upper surface 


of the slab. 


Marine gravity data produced by Hayes (1966) indicates 
that there is a small positive anomaly of approximately 20- 
30 milligais extending 200-500 km toward the continental 
side of subduction across the Chile trench at 23°S. This 


anomaly is consistent with Fig. 20 (curve D) fora slab 
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subducting at 8.0 cm/year. The subduction velocity and angle 
of descent are consistent with data available on the Chile 
trench (Le Pichon et al, 1973) which indicates the 
lithosphere plate velocity to be about 6.0 cm/year. The 
observed negative anomaly of approximately -200 milligais, 
in the Chile trench, is due to low-density ocean-floor 
sediments near the trench axis (Grow and Bowin, 1975). Watts 
and Talwani (1973) calculated the gravity effects in 
Subduction zones and concluded that the gravity effect of 
the downgoing slab is confined to the island arc and trench 
region, which is consistent with gravity profiles calculated 
here. The gravity profiles considered allow some 
investigation into the study of slab dynamics. A broad 
regional negative anomaly (Fig. 20, curve B) may indicate 
that a substantial amount of partial melt has risen from the 
top of the slab creating a low density region over the 
continental side of subduction. The anomaly can also 
indicate that the slab is moving with a relatively small 
velocity (approximately 1.0 cm/year) preventing large 
temperature gradients from occurring allowing the slab to 
warm to mantle temperatures). Alternatively, a large 
positive anomaly may indicate a slab moving at greater 
subduction velocities (approximately 8.0 cm/year) not having 


enough time to warm to mantle temperatures. 


Figure 21 gives results for different slabs dipping at 
26.6° with various heat sources. Figure 21 (curve A) shows 


the gravity effect due to a slab with shear heating and 
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Gravity profiles across the horizontal extent of 
subduction for a 26.6 dipping slab. Curves A and 
B correspond to the models of Figs. 13 and 15 res- 
pectively. Curve C is associated with a slab des- 
cending at 12.5 cm/year with only heat conduction 
considered. Curve D is associated with a slab des- 
cending with the same velocity but with adiabatic 
compression and shear heating added. 
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adiabatic compression having a_ vertical velocity of 0.7 
cm/year (aS in Fig. 13). The resulting anomaly is negative 
over the trench region, moving to near zero 400 km to the 
left of subduction (as in the figure) and falling to -40 
milligals about 1000 km from the point of subduction on the 
continent side. Figure 21 (curve B) represents the gravity 
effect of a slab moving with a vertical velocity of 0.7 
cm/year with rising melt (as in Fig. 15). There is a large 
negative anomaly with a trough of -140 milligals due to the 
low density rising material. The effect of the cold Sinking 
Slab is masked by the rising melt except for a short 
wavelength anomaly at the 1000 kn point along the horizontal 
axis. Figure 21 (curve C) is the gravity effect due to a 
Slab dipping at 26.6°, moving with a vertical velocity of 
5.6 cm/year with only heat conduction considered. Figure 21 
(curve D) is similar to Fig. 21 (curve C) except for 
adiabatic compression and shear heating considered as heat 
sources. The resulting positive anomaly peaking at 150 
mMilligals in Fig. 21 (curve C) is fovered by the shear 
heating and adiabatic compression as shown in Pig. 21 (curve 


Dh. 


Comparison of Figs. 20 and 21 indicates that the 
gravity anomalies for the shallower dipping slab are of much 
larger amplitude and extend over a larger region. This is 
due to. the descent mechanism which causes the shallower 
dipping slab to maintain a greater thickness and therefore a 


larger region of density anomaly. The implications for slab 
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dynamics are that a shallower dipping slab (ieee. 26.69) 
which has penetrated into the mantle to the same depth as a 
steeper dipping slab (i.e. 45°) will Gisplay a larger 
gravity anomaly due to the descent mechanism involved in 
Sees i as well as the velocity of subduction over the 
Same time periods considered. The descent mechanism causes a 
Slab dipping at._a shallower angle (i.e. 26.6°) to maintain a 
greater thickness in the direction of subduction are well as 
being subducted over a greater horizontal region to reach 
the same depths as a steeper dipping slab (i.e. 45°). The 
result is more slab material being subducted in the 
Shallower dipping slab adding to the gravity anomaly. The 
greater subduction velocity in the shallower dipping slab 
causes the isotherm to be pulled further into the mantle 
creating a larger positive anomaly than a_ steeper dipping 


Slab, over the same time period. 


3.5 Summary of Results 


The results presented in this chapter describe the 
thermal and gravity effects of a downgoing slab heated by 
conduction, adiabatic compression and viscous dissipation. 
From the results, the effect of dip angle on the thermal 
regime is to maintain lower temperatures with the slab at 
gentler dips (i.e. 26.69) over the same time period and 
depth interval considered as compared with a steeper dipping 
Slab (i.e. 45°). This result may be explained by the greater 


subduction velocity needed to reach the same depths in the 
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Same time periods considered as well as the greater slab 
thickness requiring more energy to be heated to surrounding 
mantle temperatures in the gentler Sloping slab. This effect 
is largely negated along the upper surface of the slab when 
Shear heating is considered. Shear heating raises the 
isotherms locally to a greater extent in the gentler sloping 
Slab due to the larger region over which this heating 
occurs. The descent of the slab causes a dip in the heat 
flow near the point where subduction begins resulting fron 
the inability of the slab to reach thermal equilibrium with 
Surrounding material at shallow depths. The various heat 
sources within the Slab do not affect the heat flow at the 


surface in the time periods considered. 


The Sib de eave velocity affects the temperature field 
within and around the descending slab. A large subduction 
velocity (i.e. 8.0 cm/year) causes the temperatures within 
the slab to remain much lower than the surrounding mantle as 
compared with a low subduction velocity (i.e. 1.0 cm/year). 
This effect, in turn, requires a greater amount of shear 
heating in the faster moving slab to create partial melt 
along the upper surface of the slab. The faster moving slab 
creates a lower minimum in the heat flow profile near the 
point where subduction begins as compared to a slower moving 


Slab. 


Rising melt, generated at depth, creates a high in the 
heat flow pattern over the horizontal extent to which this 


melt occurs. The region of rising melt is larger in the 
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Shallower dipping slab (i.e. 26.69) creating a heat flow 
high over a wider region than the steeper dipping case due 
to the depth range over which the melt is forced to rise 
(the same in both cases). Large horizontal temperature 
gradients occur at the vertical boundaries of the rising 


melt. 


The density variations within the subduction zone 
resulting from the perturbed thermal regime create anomalous 
gravity effects at the earth's surface. The cold Sinking 
Slab creates a broad positive anomaly over the horizontal 
extent of subduction resulting from the slab being denser 
than the surrounding medium. The gravity effect of the 
Sinking slab is masked by the creation of partial melt at 
the surface of the slab which rises creating a negative 
anomaly over the region of interest. The effect is more 
pronounced over the Shallower sloping slab (i.e. 26.69) 
because of the larger region over which partial melt occurs 
aS compared with the steeper dipping case (i.e. 45°). These 
gravity effects may be altered by the presence of density 


anomalies in the crust as well as topography variations. 


The variations in heat flow and gravity within a 
Subduct ion zone presented in this chapter have implications 
in the study of slab dynamics. By studying the changes in 
the observables, it may be possible to obtain a better 
understanding of processes occurring within subduction zones 
and quantitative estimates of the magnitudes of these 


processes. An attempt has been made to determine changes in 
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heat flow and gravity as a result of a slab descending into 
the mantle under various conditions. It is hoped that the 
results obtained may be applied to the further study of slab 
dynamics in regard to heat flow and gravity in subduction 


zones and related island arc regions. 
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CHAPTER 4 THE EFFECT OF SUBSURFACE TEMPERATURE 
VARIATIONS ON THE BEHAVIOUR OF TIME-VARYING 


ELECTROMAGNETIC FIELDS 


The electrical conductivity of the earth is in general 
a complicated function of several variables including 
temperature and pressure. The effect of pressure on 
conductivity is not well known although some work has’ been 
done on the variation of conductivity with pressure-induced 
phase transitions at depths between 400 km and 900 km 
(Rikitake, 1959, Akimoto and Fujisawa, 1965). Although the 
present numerical analysis includes regions at depths of 600 
km, the variation in conductivity as a function of pressure 
is ignored. The relationship between temperature and 
electrical conductivity can be expressed as: 


-E,/2kT -E,/2kT -E,/kT 
Ovo, "ve + To" e + Gi? eo i (14) 


taking into consideration three different conduction 
mechanisms; impurity, intrinsic, and ionic (Rikitake, 1966). 


In the above equation, (o"), denote the conductivities 
9 G 


293 

at infinite temperature, (E) , ht denote the excitation 
> > 

energies, T is the absolute temperature and k is the 

Boltzmann constant. The temperature dependence of the 


excitation energy is contained in the values for (E); . This 


dependence occurs because of the dominance of ditterent 
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conduction terms in specific temperature intervals. Although 
the excitation energy is expected to depend upon pressure, 
only crude estimates have been obtained as to an exact 
relationship (Rikitake, 1966) and therefore the effect has 
been ignored in the present models. Since the composition of 
the mantle consists mostly of olivine-pyroxene minerals, it 
is reasonable to choose the following values for the 
constants: 0, = 107%,o2" = 103,03! = 10-6, Pia bo tao 6 
where of is in emu and E is in electron volts (Rikitake, 
1966). The conductivity distribution does not take into 


consideration specific variations related to melted material 


other than those due to temperature. 


ia 


The thermal models considered in this work are: 


Model (a): This model gives the temperature distribution 
corresponding to the time when the slab just begins to 


descend (Initial state), 


Model (b): This gives the temperature distribution after 
the slab has descended for a period of 21.2 million years 
with heat conduction considered as the only mechanism of 


heat transport (Final state), 


Model (c): The period of subduction is the same as that of 
model («b), but the temperature distribution is influenced by 


partial melting along the upper surface of the slab, 
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Model (d): The period of subduction is the same as that of 
model  (b), but this model assumes that the material melted 
along the surface of the slab rises to within 30 km of the 
earth's surface. In models (b), (c), and (d) the angle of 
subduction of the slab is 45° and its downgoing velocity is 
1.0 cm/year. The vertical thickness of the slab is taken to 


be 80 km. 


Figure 22 shows the temperature distribution obtained 
for model (a) (Fig. 22a - left) and for model (d) (Fig. 22b 
- left). The temperature distributions for models (b) and 
(Cc) are not shown. The dimensions of each model are 300 km x 
300 km. The temperature varies from 0°C at the surface to 
approximately 1700°C at a depth of 300 km. The electrical 
conductivity model shown in Fig. 22a (right) corresponds to 
the thermal model (a), and that shown in Fig. 22b (right) 
corresponds to model (d). Since the finite-difference 
numerical method for the electromagnetic calculations 
requires that the conductivity distribution be divided into 
various cells, each with a specified conductivity value, 
small variations in temperature (and conductivity) are 
ignored. The values assigned to each conducting cell are 


Shown in Fig. 22. 


Figure 23 shows the geophysical model considered in 
this work. It consists of a conducting layer of thickness 
‘d* and of two different conductivities o, and o, overlying 
the conductivity distribution determined by each of the four 


thermal models. The inducing source field is assumed to be 
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Temperature distribution and electrical 
conductivity models for (a) the initial 
state, and’ (bb) the final state. .The 
conductivity values are in emu units. 
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The geophysical model. 


7 


Bowe: s 


~ ad 


, 2 
ewe 
ee 
o 36 


ow 
galt = 

i ried 

4 4 ay 
7a oe 


; anes ‘ 
et sate os 
oe Ok 


qe es A 
othe Ste 8g 
b8? 50° Big 


) 


y . , y . : A e § iG Ff fi . i ; | | : a 
‘9, Pehom Isotoumigoga oat VER sgt 9 8 eke sees 


P ‘ ; 
ee re mt om" Te) 
} 1 oy = i, H aati eek? BRED *) = 4 aut 


70 


uniform with its horizontal electric component in the x- 
direction (E-polarization). The frequency of oscillation is 


chosen to be 0.75 Hz. 


The electromagnetic modelling technique used in the 
present work is basically the same as that described by 
Jones and Pascoe (1971) and Pascoe and Jones (1972), with 
the modifications. proposed by Williamson et al. (1974) and 
Brewett-Taylor and Weaver (1976). The technique is based on 
the finite-difference approach to the solution of Maxwell's 
equations in two aidenstens and has been extensively used in 
investigating various conductivity anomalies (Jones, 1973). 
In the numerical method a 71 x 71 grid of mesh points with 
variable grid dimensions is superimposed on the geophysical 
model shown in Fig. 23. It is noted that, for large positive 
and negative values of y, the conductivity models of Figs. 
22 and 23 have layered structure, which is one of the 
boundary conditicns to be satisfied in the numerical method 
of Pascoe and Jones (1972). It should be mentioned here that 
for the convenience of display only a part of the 71 x 71 
grid used is shown in Fig. 22. The lower boundary in the 
electromagnetic mumerical model has been chosen at z = 630 
km, whereas the lower boundary in the heat flow naumerical 
model is at z = 300 km. This difference in the depths of the 
lower boundary implies that it has been assumed the 
temperature is constant between depths of 300 km and 630 kn. 
This assumption permits the lower boundary in the electrical 


conductivity configuration to be far away from the earth-air 
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interface z = 0. 


4.2 Discussion of Results 


a —S 


The numerical calculations for the electric and 
magnetic field components have been carried out for two 
different cases: (1) the conductivity of the lower layer is 
varied while the conductivities of the upper 
layer, 91 and o2, and its thickness are kept constant, and 
(2) the conductivity o2 is varied while keeping d, o, and 
the conductivity of the lower layer constant. In the 
following discussion, E, denotes the amplitude of the 
horizontal electric field component, whereas Hy and Hy, 
denote the amplitudes of the horizontal and vertical 
Magnetic components, respectively. The apparent 
resistivity is defined as ois (0.2/f) (Ex/Hy)?2, where f is 
the frequency in hertz, Ey is in nV/ka, Hy is in oersteds 
and P. is in emu. The surface heat flow is expressed in 
units of microcalories/ (cm@-sec). To distinguish between the 
electromagnetic models and the heat flow models, letters 
A,B,C and D refer to the electromagnetic models (Fig. 23) in 
which the lower layer configuration is determined by the 
conductivity distribution corresponding to the heat flow 


models a,b,c or d, respectively. 


To study the effect of the subsurface temperature 
variations on the electric and magnetic fields observed at 
the earth's surface, calculations were carried out for 


models A,B,C and D by choosing d= 5 km for the upper layer 
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Fig. 24. The spatial variation of heat flow and the 
electromagnetic field components at the 
surface z=0 for different conductivity 
models, with d=5 km and 0,=02=10 !° emu. 
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depth and 91 = 602 = 10-45 emu for the upper conductivities. 
Figure 24 shows the spatial behaviour of the amplitudes E,, 
Hy, Hz and their ratios for points at the surface z = 0 
along the y-direction (i.e. perpendicular to the source 
electric field). Although the results are presented only for 
the region -150 km < y* < 150 km, the left and right outer 
boundaries in the © numerical calculations for the 
conductivity models are kept far beyond y = + 150 km. This 
is done in order to satisfy the boundary conditions which 
require that in regions’ far away from the anomaly, the 
conductivity structure is layered and that the fields there 
remain unperturbed. The factors that appear on the ordinates 
of some of the amplitude curves are the scaling factors 
which are to be used in determining the exact amplitude. For 
example, at y = -150 km, the apparent resistivity for 


model A has the value of approximately 1.0x1015 emu's. 


It is evident from the results shown in Fig. 24 that 
the surface heat flow F for the final state (models B-D) is 
less than that for the initial state (model A). This is 
mainly due to cooling in the mantle since the beginning of 
subduction. AS the high temperature isotherms rise to the 
surface (models C and D) F shows a general enhancement. The 
heat flow F, which undergoes an abrupt change near y = 40 kum 
for the initial state, exhibits a smooth variation for the 
final state, with a slight horizontal shift of about 5 km in 
the location of itS minimum. The ‘kink that appears in F 


for the model D shows the effect of rising hot material from 
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the upper surface of the descending slab. Since the results 
for the electromagnetic field components shown in Fig. 24 
correspond to the electromagnetic model with uniform source 
field and a uniformly conducting overburden, any variation 
in the spatial behaviour of the field components can be 
attributed to the subsurface temperature variations. The 
results for the horizontal electric component indicate that 
Ex in general has low amplitude in the regions where there 
is increased heat flow and high amplitude in the regions 
with decreased heat flow. Noting that the only difference 
between models B and C is that the latter assumes the 
presence of melted material along the slab, it is readily 
observed that these materials which are responsible for the 
increased heat flow cause a reduction in the Ex amplitudes 
observed at the surface. The rising of hot material from the 
Slab (model D) considerably affects both the amplitude and 


Shape of the Ey, component. 


AS expected, the horizontal magnetic component Hy has 
large amplitude in regions of increased heat flow. In the 
case of the initial temperature distribution (model A), Hy 
increases gradually as the region where F undergoes a sudden 
variation (near y = 40 km) is approached, reaching a maximum 
near y = -10 km over the higher heat flow side. Beyond this 
point, Hy falls off and gradually approaches a constant 
value for large values of j{y{ aS required by the uniform 
source field. For model B however, Hy does not exhibit any 


pronounced maximum but instead undergoes a hininum 
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approximately where the surface heat flow is a minimum. As 
the high temperature isotherms rise to the surface (models C 
and D) Hy again Shows behaviour similar to that of Hy for 


model A. 


The behaviour of the vertical magnetic component (H,) 
at the surface eee to be more sensitive to the 
subsurface temperature variations than the behaviour of the 
horizontal. component (Hy)- This is particularly evident in 
the horizontal shift in the location of the peak value of Hz 
which is considerably affected by the rate of decrease (as a 
function of y) of the surface heat flow. A comparison of the 
behaviour of Hz with that of Hy indicates that with respect 
to the location where F undergoes a rapid variation, the 
maximum amplitude for Hy always occurs over the higher heat 
flow regions whereas the maximum of Hz occurs some distance 
to the right over the decreased heat flow side. The results 
also indicate that the sovement of hot magma from the _ slab 
(model D) has considerable influence on the amplitude as 


well as on the general behaviour of the vertical component. 


The results for the vertical to horizontal magnetic 
field ratio (Hz/Hy) and the apparent resistivity (p.) are 
also given in Fig. 24. The behaviour of Hz/Hy as a_ function 
of y is very Similar to that of HZ. The apparent resistivity 
is, as expected, less over the higher heat flow regions than 
over the lower heat flow regions. Like the other field 
components, P, is also influenced by the subsurface 


temperature variations, the most noticeable effects being 
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due to the melted material that has risen from the Slab. 


Figure 25 shows the contours of equal amplitude of the 
subsurface temperature distribution and of the three field 
components for models A-D. Only the temperature distribution 
below the surface z=0 is shown, since it has been assumed 
that in the numerical model for the heat flow problem, the 
temperatures for z < 0 remain unperturbed and are taken to 
be zero for convenience. For the electromagnetic fields, the 
amplitude contours include both the free-space region 
(zZ < 0) and the conducting region (z > 0). It must be noted 
that the contour values of each of the three components are 
the same for the four models. One feature that is readily 
apparent is that the effect of the subsurface temperature 
variations is more pronounced in the amplitude contours of 
Hy and Hz than in the contours of Ey. The contours for the 
vertical component for the four models show clearly the 
horizontal shift in the location of the Hz maxima that was 
Observed earlier (Fig. 24). The presence of the two maxima 


in Hz, is also evident in the contour plots for model D. 


The effect on the electromagnetic fields of varying the 
upper layer conductivity contrast has been studied by 
varying the conductivity o2while keeping both 0:1 and d 
constant. Figure 26 gives the results from H,/Hy and p, at 
the surface for various values of 62/0, for model D 
with d= 2 km and 0)= 107-45 emu. At the frequency of 0.75 
Hz, the skin depths in the conductivity o2 are 58, 18, 5.8, 


and 1.8 km for 02/0; = 1,10,102, and 103, respectively. The 
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Fig. 26. The spatial variation of Hz/Hy and (0, at the 


surface z 


QO for different upper layer 


conductivity contrasts for model D, with 
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and 02 is located at y=47.5 km.) 
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boundary between 0); and J2is located at y = 47.5 km. For 
convenience of display, different scale factors are used in 
the results for Hz/Hy. For example, taking into account the 
scale factors shown in Fig. 26, the maximum amplitudes of 
Hz/Hy should read as Ustad, Onegn (1.5 and “2.0 for 
O2v0O1 = 1,10,102, and 103 respectively. The behaviour of 
Hz/Hy for the model in which o2/o, = 1 is similar to that 
Shown in Fig. 24 for model D. Since, in this case, the upper 
layer is unifcrmly conducting, the variation in the field 
components as a function of y can be interpreted as due to 
Subsurface temperature variations only. The results 
indicate, that as 02 increases, the effect of the upper 
layer conductivity contrast increases and the effect of 
temperature variations is dominated by the upper layer 
contrasts. This is evident from the sudden increase in the 
amplitude of Hz/Hy near y = 47.5 km for large values of o2. 
Since 02/0, = 103 is of the same order as that for the sea- 
continent interface, the results for 02/o, = 103 would be 
applicable to the sea-land model in which o2 represents the 
seae Thus, from these results it is apparent that large 
conductivity contrast at the surface dominates the behaviour 
of the magnetic field ratios and the apparent resistivity to 
the extent that the effect of the subsurface temperature 
variations is suppressed. Although the effect of the surface 
conductivity contrast will in general be different for 
different frequencies, its predominance in this modei may be 
due to its location (y = 47.5 km) which is not far from the 


region qnear y = 35 km) where the surface heat flow 
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undergoes a rapid variation. 


4.3 Summary of Results 


The results presented in the previous section describe 
the behaviour of time-varying electromagnetic fields 
corresponding to heat flow distributions associated with a 
downgoing slab. From the results presented for a uniformly 
conducting upper layer several aspects of slab dynamics may 
be investigated. In the case of a descending slab heated 
only by conduction, the results indicate considerable 
variation in the vertical magnetic component (Hz) as well as 
the magnetic field ratio Hz/Hy. The effect of the cold 
descending slab is to increase both of these quantities over 
the region of slab descent. The effect of shear heating 
along the upper surface of the slab causes a slight increase 
in the value of the horizontal magnetic component Hy with 
little effect on the other field components. The inclusion 
of rising melt from the upper surface of the slab causes 
Significant variations in the electric and magnetic field 
components with the largest variations observed in the 


vertical magnetic component (Hz) and the magnetic field 


ratio Hz/Hy- 


From the above results the study of electric and 
‘magnetic fields in subduction zones may be used to. study 
Slab motion. A decrease in the apparent resistivity may be 
indicative of melted material rising from the top of ‘the 


descending slab as well as corresponding increases in the 
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magnetic field components. nite effect of a slab descending 
with no melting along the upper surface is to increase 
apparent resistivity while decreasing the horizontal 
Magnetic component. When a lateral conductivity contrast is 
introduced in the upper layer, the results show that the 
effect of the temperature variation in the lower layer is 
reduced, and that, for a large conductivity contrast (e.g., 
sea-land interface), the subsurface temperature variation 
has little effect on the surface electromagnetic fields. 
Although this result is expected on the basis that currents 
induced in good conductors flow near the surface, in this 
model, the dominance of the surface conductivity contrast 
over the subsurface temperature variations is probably due 
to the closeness in the location where the heat flow varies 
rapidly. This aspect of the screening effect of the surface 
conductivity contrast as well as the effects on the 
electromagnetic fields of different frequencies requires 


further investigation. 
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